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Radiative  Transfer  in  an  Atmosphere-Ocean  System: 

A Matrix  Operator  Approach 

George  W.  Kattawar,  Terry  J.  Humphreys,  and  Gilbert  N.  Plass 

ABSTRACT 

It  is  the  purpose  of  this  paper  to  demonstrate  how  the  matrix  operator 
method  can  be  effectively  implemented  to  couple  the  radiation  fields  of  the 
atmophere  and  ocean.  Azimuthal ly  averaged  radiances  and  irradiances  are  presented  as 
a function  of  optical  depth  for  a conservative  Rayleigh  scattering  medium  of  total 
optical  thickness  = 1000  with  a dielectric  interface  placed  at  optical 
depths  of  0.01,  0.1,  1.0,  and  10.0,  and  for  various  solar  incident  angles. 

A INTRODUCTION 

Mo- v 

In  order  for  us  to  understand  the  effect  of  changes  in  optical  properties 
of  either  the  atmosphere  or  ocean  on  the  radiation  field,  we  must  be  able  to 
perform  accurate  radiative  transfer  calculations  for  this  coupled  system.  It 
is  important  for  researchers  in  the  remote  sensing  area  to  know  not  only  what 
regions  of  the  spectrum  are  affected  but  to  what  degree  they  are  affected  by 
known  changes  in  the  optical  properties  of  the  medium.  One  area  of  extreme 
importance  is  the  remote  sensing  of  phytoplankton  in  a body  of  water  since 
these  organisms  essentially  determine  the  effective  productivity.  In  addition 
to  the  radiation  in  vitro  for  remote  sensing  purposes  one  is  also  interested 
in  the  radiation  in  situ  for  biological  purposes  since  the  irradiance  is  an 
important  factor  in  photosynthesis. 

To  perform  a radiative  transfer  calculation  at  any  particular  frequency, 
we  need  only  two  quantities,  namely  e(o)  (the  volume  scattering  function 
in  units  of  m’1  sr"1)  and  c (the  attenuation  coefficient  in  units  of  m'1). 

Now  if  the  medium  is  inhomogeneous,  then  we  must  specify  these  quantities 
at  each  depth  within  the  medium.  Other  useful  quantities  derivable  from 
these  two  are: 
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b (scattering  coefficient)  = 2 ^/^(o  ) sin 


P (o)  (phase  function)  = B(o)/b, 


where 


2uo/n  P( o ) sin  0 do  = 1 , 


a (absorption  coefficient)  = c-b 


L:  (single  scattering  albedo)  - b/c , 


i (optical  depth  between  layers  at  z,  and  z~)  - /.  c dz.  (61 

i z z] 

It  should  be  emphasized  that  all  of  the  optical  parameters  as  given  are 
for  all  constituents  within  the  sample  being  considered. 

If  we  let  I(x,u,<t)  denote  the  radiance  then  the  equation  of  transfer 
assumes  the  following  form 


(p  +1)1  (i  ,p»<J>)  = /|.|  P(t,u,<|>;p'  ) I (t  ;p;4>'  )du'd<t 1 (7) 


where  we  have  adopted  the  convention  that  p > 0 denotes  downward  moving  radia- 
tion and  v < 0 denotes  upward  moving  radiation.  It  is  the  purpose  of  this 
article  to  demonstrate  that  quite  accurate  solutions  to  this  equation  can  be 
obtained  by  the  matrix  operator  technique  subject  to  various  boundary  condi- 


tions. 


BASIC  THEORY 


The  historical  development  of  the  matrix  operator  theory  can  be  found 
in  Plass,  et.al.1,  and  we  will  not  go  into  the  details  of  discretization 
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or  Fourier  decomposition.  The  fundamental  problem  arises  when  we  couple  the 
atmosphere  to  the  ocean  through  a dielectric  interface.  Let  us  assume  that 
we  have  a discrete  set  of  N cosines  of  the  polar  angle  0 < 0 < tt/2  for  the 
ocean  such  that  = cos  o.  > 0.  Similarly  we  will  have  a set  of  points 
VK  = -cos  < 0 which  will  label  the  upper  hemisphere  (see  Fig.  1).  Out 
of  the  set  of  N points  for  the  ocean  only  M of  them  can  be  mapped  into  the 
atmosphere  i.e.,  those  for  which  p.  > cos  e where  6 is  the  critical  angle 

I v*  v 

defined  by  sin  0c  = 1/n  where  n is  the  refractive  index  of  water.  The 
remaining  N-M  points  arc  restricted  to  the  ocean  and  cannot  be  mapped  into 
the  atmosphere.  These  points  specify  the  region  of  total  internal  reflection. 
The  relationship  connecting  « to  o is  simply  Snell's  law,  namely 


or 


sin  a..  = n sin  (m 


(8) 


^ = (1-  0-u*)n2)1/2 


(9) 


where  c.  - cos  . The  selection  of  the  quadrature  mapping  is  the  most 
crucial  aspect  of  the  matrix  operator  approach.  The  constraints  placed  on 
the  selection  are  as  follows.  First  we  must  be  able  to  properly  normalize 
a phase  function  for  the  atmosphere  due  to  the  aerosols  as  in  Eq.  (3). 
Secondly  we  must  be  able  to  conserve  the  energy  reflected  from  and  trans- 
mitted through  the  interface.  Finally  we  must  be  able  to  properly  normalize 
a highly  asymmetric  phase  function  due  to  the  hydrosol  scattering  in  the 
ocean.  Our  solution  to  this  problem  has  been  enormously  successful  and  it 
proceeds  as  follows.  We  first  consider  a Gaussian  quadrature  of  order  M 
for  the  atmosphere  which  is  mapped  into  the  half-space  i.e.,  0 < ? s 1 
and  of  course  0 < <j>  < 2u.  The  order  of  the  quadrature  used  Is  determined  by 
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the  condition  that  it  must  be  sufficient  to  normalize  the  aerosol  phase 
function  to  within  0.1%.  To  handle  the  second  condition,  consider  a 
radiance  stream  I (0,4.)  striking  the  interface  at  some  angle  a and  entering 
the  ocean  at  the  refracted  angle  e.  For  the  case  of  a continuum  of  radiation 
1(0,4-)  striking  the  interface  from  above, the  radiance  just  below  the  interface 
is  given  by 

1(0,4.)  = !(«,*)  t (a)  (10) 


where  t{«)  = 1 - r(u)  and  r(u)  is  the  Fresnel  reflection  coefficient  for 
unpolari/.ed  radiation  given  by 


r(o)  = [tan^(a-o)/tan^(o+n)  + sin^(o-o)/sin^(u+o)]/2 


(11) 


This  result  holds  for  oblique  incidence,  i.e.  a f 0.  For  normal  incidence 
Eq.  (11)  reduces  to 

r(0)  - ("  i j)2  (12) 

Using  Eq.  (8)  , Eq.  (10)  can  be  written  as 


I(o, 4-)  = I(u,4>)  t (a)  n^  (13) 

2 

a result  first  noted  by  Gershun  . Mow  in  the  discrete  approximation  Eq.  (10) 
becomes 


l(o., t) 


= I(o.  ,$)t(«.) 


Cf 


1 ’■  "t  c? 


(14) 
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A 0 

where  and  are  the  weights  for  the  atmosphere  and  ocean  respectively 
associated  with  the  quadrature  integration.  Comparing  Eq.  (14)  with  Eq.  (13) 
we  can  now  determine  the  ocean  weights  C?  uniquely,  namely 


c° 

i 2 

V 


(15) 


2 2 1/2 

Now  ni  = (1  - (1  - S^)/n  ) where  0 < < 1 and  pc  < p..  < 1 where 

2 1/2 

pc  = (1  - 1/n  ) . The  function  c/p  is  shown  in  Fig.  2b.  The  important 

aspect  of  this  function  is  that  it  is  bounded  and  has  bounded  derivatives 
of  all  orders.  It  is  ideally  suited  for  Gaussian  quadrature  integration 

since  the  error  term  for  an  nth  order  Gaussian  quadrature  is  proportional 
to  the  2nth  derivative  of  the  function.  Since  the  phase  function  is  decom- 
posed in  a series  of  Legendre  polynomials,  we  must  be  able  to  integrate  these 

polynomials  with  high  accuracy  wit!)  our  mapped  points  and  weights  for  the 

1 k 

ocean.  Let  us  now  consider  f p dp  for  the  ocean:  where  k is  an  integer. 

c 

In  the  discrete  approximation  this  becomes 


M k r0  1 
c'  „ 


v k-]  r rA 

2 ih  wf  h ci 


(16) 


It  should  be  noted  that  if  k is  odd  then  Eq.  (16)  shows  that  we  ere  effectively 

integrating  a polynomial  of  degree  k in  the  atmosphere  and  this  result  will  be 

exact  as  long  as  k < 2m-l . However  if  k is  even  then  we  no  longer  have  a 

finite  polynomial  since  the  function  we  are  integrating  is  proportional  to 
? 2 ( k-1 ) 12 

C[1  - (1  - c )/n  j . However  this  function  is  quite  well  behaved  and  has 

bounded  derivatives  of  all  orders  and  again  Gaussian  quadrature  should  do 
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quite  well.  We  have  tested  these  ideas  using  an  8th  order  Gaussian  quadrature 


and  for  0 < k < 15  the  worst  error  we  encountered,  which  was  for  k = 0,  was 
1 2 

2 parts  in  10  . All  odd  values  of  k were  exact.  Even  for  k = 49  the  error 

was  only  2 parts  in  104.  The  next  logical  question  to  be  asked  is  what  errors 

are  encountered  if  one  uses  Gaussian  points  within  the  critical  angle  for  the 

ocean  and  maps  them  into  the  atmosphere?  This  was  the  technique  used  by 
3 

Tanaka  and  Nakajima  . We  will  now  show  why  this  method  has  serious  defects. 

A 

With  this  method  the  atmospheric  weights  are  determined  from  the  ocean 
weights  C*?  by 


CA  = V 


C° 
' 1 


(17) 


where  the  ocean  weights  are  determined  by  using  Gaussian  quadrature  mapped  into 

1 k 

the  region  pc  < p < 1 . If  we  again  consider  /q  f ; d{,  then  in  the  discrete 
approximation  we  got 


M 

>: 

i~l 


M 

y 

i-1 


k-1  r0 
l . n . L . 
'i 


(18) 


For  k = 0 we  show  a plot  of  the  function  p/f,  as  a function  of  p in  Fig.  2a. 

Not  only  is  the  function  unbounded  at  p = p but  it  lias  unbounded  derivatives 
at  this  point  and  therefore  quadrature  integration  should  produce  large  errors. 
In  general  one  will  be  integrating  p[l  - n (1  - p )]'  ' which  should  be  “act 

for  odd  values  of  k up  to  2m-l . However  for  k=2  the  function  becomes  bounded 
but  it  still  has  unbounded  derivatives  and  quadrature  will  produce  relatively 
large  errors.  For  even  values  of  k > 2 the  function  is  bounded  and  has  bounded 
derivatives.  We  tested  this  scheme  with  an  8th  order  Gaussian  quadrature. 

5 

For  k = 0 the  error  was  4.6%  and  for  k = 2 it  was  1 part  in  10  and  becomes 
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better  for  higher  values  of  k(even).  The  error  for  k = 49  was  of  the  same 
order  as  the  earlier  method  described,  namely  2 parts  in  104.  It  should  now 
be  clear  that  this  latter  method  would  be  particularly  bad  if  one  is  considering 
Rayleigh  scattering  for  the  atmosphere  which  is  the  dominant  mode  of  scattering 
for  the  earth's  atmosphere  in  the  visual  region  of  the  spectrum. 

MATRIX  OPERATOR  CONSTRUCTION 

We  will  now  show  how  to  construct  the  appropriate  matrices  from  the  opera- 
tors associated  with  the  atmosphere,  ocean,  and  interface.  Throughout  the 
sequel  we  will  only  consider  azimuthally  averaged  quantities  which  involve  only 
the  first  term  in  the  Fourier  decomposition  of  the  phase  function  (see  Plass, 

et.alJ)  . The  extension  to  the  complete  decomposition  is  straightforward . 

A A A 

Let  k U-.C-)  = R..  = R be  the  matrix  representation  of  the  reflection  operator 
1 J 1 J 

R^  for  the  atmosphere  for  radiation  incident  on  top  of  the  layer.  This  matrix 
has  dimension  MxM  where  M is  the  order  of  quadrature  for  the  atmosphere.  Simi- 
larly is  the  matrix  associated  with  the  transmission  operator  for  the  layer. 
It  also  has  dimensions  MxM.  for  the  ocean  without  the  interface  we  have  the 
corresponding  matrices  R^(  w i , w ^ ) = r9.  = R^  and  T®(ij.,ji.)  = T?.  = T°  of 
dimensions  NxN.  The  reflection  matrix  R^w  going  from  air  to  water  is  diagonal 
and  simply  consists  of  the  Fresnel  reflection  coefficients  defined  by  Eq.  (10). 
Similarly  the  transmission  matrix  from  air  to  water  T^  = E - R^w  where  E is 
the  identity  matrix  and  both  have  dimensions  MxM.  However  to  formalize  our 
matrix  products,  T^w  must  have  dimensions  NxM,  we  must  insert  zeros  for  all 
elements  in  rows  M+  1 through  N.  Now  the  reflection  matrix  from  water  to 
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air  Rw^  is  again  diagonal  and  the  first  M elements  along  the*  diagonal  are 
again  the  Fresnel  reflection  coefficients  whereas  the  remaining  N-M  elements 
are  all  unity  corresponding  to  total  internal  reflection.  Similarly  the 
transmission  matrix  from  water  to  air  is  TWA  = E - R^.  It  should  be  noted 
that  this  matrix  has  all  zero  entries  in  rows  M+I  through  N.  Therefore 
this  matrix  can  be  reduced  to  dimension  MxN.  We  can  now  use  the  matrix 
operator  formalism  to  set  up  the  reflection  and  transmission  matrices  for  the 
combined  system,  namely,  atmosphere,  interface,  and  ocean. 

To  illustrate  the  adding  algorithm  (see  Plass1,  et.al.)  the  matrices  for 
the  combined  interface  and  ocean  are 


-1 


-10  = -AW  + W-  ' 'OW  -O'. 


AW 


*10 


= V~  ‘ -WA'O^ 


*AW 


(19) 

(20) 


It  should  be  noted  that  the  matrix  operations  encountered  in  Eqs.  (19)  and  (20) 
are  at  least  dimensionally  consistent.  The  matrix  Rjq  has  dimensions  MxM 

whereas  the  matrix  Tjq  has  dimensions  NxM.  The  physical  interpretation  of 
equations  (19)  and  (20)  is  extremely  sin.pl c if  wc  consider  writing  the  expression 
(E-RqRw^)  1 formally  as  a power  scries,  namely 


,-l 


~ + -nBu/i  + Bn^uA-n-WA  + • • • 


0 WA  '0  WA  ' O.'WA 


(21) 


Therefore,  using  the  first  term  in  Eq.  (21)  and  reading  Fq.  (19)  from  right 
to  left,  we  see  that  T^  transmits  the  radiation  through  the  interface,  next 
Rq  reflects  it  from  the  ocean  and  finally  T^  transmits  it  back  through  the 
interface  and  into  the  atmosphere.  Now  the  second  and  higher  order  terms  in 


Eq.  (21)  simply  give  all  orders  of  reflection  between  interface  and  ocean 
before  emergence.  Eq.  (20)  has  a similar  interpretation.  It  should  also  be 
pointed  out  that  one  can  also  obtain  the  radiation  field  at  any  interior 
point  in  the  medium  simply  by  having  the  reflection  and  transmission  matrices 
for  the  medium  both  above  and  below  the  interior  point. 


PHASE  FUNCTION  TRUNCATION 


One  of  the  major  difficulties  encountered  when  solving  the  transfer 
equation  in  the  ocean  is  the  strong  asymmetry  in  the  phase  function  from  the 
hydrosol  scattering.  To  perform  a normalization  on  such  a strongly  peaked 
function  would  require  quadratures  of  extremely  high  order.  One  therefore 
has  to  consider  the  possibility  of  sacrificing  some  accuracy  for  savings 

4 

in  computer  costs.  We  have  adopted  a technique  used  by  Potter  to  eleviate 
some  of  the  difficulties  associated  with  a strongly  peaked  phase  function. 

The  technique  is  basically  the  following.  Let  P (o)  be  a strongly  peaked 
phase  function  subject  to  the  normalization  condition  in  Eq.  (3).  Also  let 
y = cos  0 . Now  in  the  delta  function  approximation  the  diffraction  peak  is 
treated  as  a Di>'ac  delta  function  and  one  approximates  P(y)  by 

P(>)  = AA(1  - y)  + P*(  ,)  (22) 

where  P'(y)  is  a phase  function  whose  diffraction  peak  has  been  removed. 

Using  Eq.  (22)  and  Eq.  (3)  we  can  evaluate  A by 

A = /]  [P(y)  - P*  (y)]dY  = - /]  P'(y)dY  . (23) 

Therefore 
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/]  P'(Y)dY  - ^ - A . (24) 

We  can  now  define  a properly  normalized  scaled  phase  function  P<.(y)  by 

PS^>  = 1-^1  A (25) 

To  insure  equality  in  the  volume  scattering  function  p.(y)  it  is  easy  to  see 
that 

bs  = b(l  - 2nA)  (26) 

where  the  subscript  S is  used  to  denote  parameters  involving  the  normalized 
truncated  phase  function.  The  optical  depth,  t,  and  single  scattering  albedo, 
wo>  transform  as  follows: 

TS  = <bs+  a)T/(b  + a)  (27) 

woS  = wo^  ' 2ttA)/ ( 1 - 2ttAwq)  (28) 

4 

This  method  has  been  tested  by  Potter  and  he  found  that  for  solar  incidence 
angles  <84°  the  irradiances  were  accurate  to  within  IX.  A comparison  between 
azimuthally  averaged  radiances  for  the  reflected  radiation  agreed  to  within 
Vi  for  angles  <87°.  The  transmitted  radiance  was  also  accurate  to  the 
same  order  except  in  the  region  of  the  incident  beam  and  for  small  optical 
thicknesses.  Once  the  optical  thickness  is  large  enough  for  sufficient 
photon  diffusion,  the  larger  discrepancies  disappear. 


g^A. 
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COMPUTATIONAL  RESULTS 

To  demonstrate  the  power  and  versatility  of  the  matrix  operator 
method  we  have  considered  a model  whose  total  optical  thickness  x , = 1000. 
We  also  assume  conservative  (wq  = 1.0),  Rayleigh  scattering  throughout  the 
region.  To  explore  the  effects  due  only  to  an  interface  we  have  placed 
the  interface  at  various  optical  thicknesses  within  this  region.  The 
interface  is  assumed  to  have  a refractive  index  of  1.338  characteristic 
of  ocean  water  at  visual  wavelengths.  The  following  optical  quantities 
are  computed  at  many  points  both  above  and  below  the  interface  and  we  will 
use  the  subscripts  U and  D to  denote  upwelling  and  downwelling  radiation 
respectively. 


I°(u)  (azimuthal  ly  averaged  radiance)  = ^ I(u,<f>)d<t>  (29) 

Hjj  Q(upward  or  downward  irradiance)  = /J  Iy  D( u ) ydjj  (30) 

hy  D(upward  or  downward  scalar  irradiance)  = jr  Iy  D(u)du  (31) 

There  are  other  quantities  derivable  from  these  which  have  important 
physical  meaning.  The  quantity  d(HD-Hy)/dx  = (1  - uQ) (hQ + hy)  is  propor- 
tional to  the  heating  rate  in  the  medium.  Also  if  the  medium  is  conservative 
( u>0  = 1.0)  then  Hy  - Hy  is  a constant  at  all  depths  and  this  is  a good  test 
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I 

of  numerical  stability  for  an  algorithm  used  to  solve  the  equation  of 
transfer.  The  quantity  hy+hy  is  proportional  to  the  energy  density  in 
the  medium.  The  ratio  of  hy/Hy  or  hy/Hy  gives  an  indication  of  the  degree 
of  anistropy  of  the  radiation  field  and  is  always  greater  than  or  equal 
to  unity.  This  result  is  obvious  from  Eqs.  (30)  and  (31).  If  the  radia- 
tion field  is  isotropic  then  hy/Hy  or  hy/Hy  = 2.0.  If  the  field  is  highly 
peaked  in  a certain  direction  say  y'  then  hy/Hy  or  hy/Hy  ~ l/u‘. 

In  Table  1 we  present  various  optical  quantities  as  a function  of 
optical  depth  for  a medium  of  optical  thickness  xmflx  = 1000  with  con- 
servative Rayleigh  scattering  throughout.  No  interface  exists  in  this 
medium.  Three  solar  incident  angles  are  considered  namely  £,q  = 0.980, 

0.592,  and  0.0199  corresponding  respectively  to  solar  zenith  angles  of 
a = 11.5°,  53.7°,  88.9°.  The  first  thing  to  note  is  that  Hy  - Hy  is  constant 
for  each  CQ  and  for  all  optical  depths  and  represents  the  downward  ir- 
radiance  exiting  the  bottom  of  the  medium.  We  have  also  assumed  a perfectly 
absorbing  bottom.  An  extremely  important  aspect  to  these  calculations 
can  be  observed  in  the  irradiances  Hy  and  Hy.  In  the  vicinity  x = 10.0 
for  = 0.980  both  Hy  and  Hy  exceed  their  values  at  x=0  by  more  than  20°', 
however  for  £0  = 0.592  a very  weak  maximum  is  reached  in  the  vicinity  of 
x = 0.1  whereas  for  = 0.0199  we  found  no  maximum.  This  irradiance 
increase  is  a purely  three  dimensional  effect  and  does  not  occur  in  one 
dimension.  A detailed  study  of  this  phenomenon  will  appear  in  a future 
publication.  The  ratio  hy/Hy  at  x=0  is  simply  1/SQ,  since  the  input  is 
monodi recti onal , but  by  x=10  sufficient  diffusion  has  taken  place  so  that 
the  internal  radiation  field  is  quite  isotropic  and  therefore  hy/Hy  = 2.0. 
The  quantity  hy/Hy  at  x = 0 shows  significant  departure  from  an  isotropic 
distribution,  however,  at  x=  10  it  shows  the  same  isotropic  behavior  as 
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the  downward  radiance.  The  reason  for  this  is  that  at  the  top  of  the 
layer  there  is  always  a significant  component  of  the  radiance  due  to 
single  scattering  which  of  course  for  Rayleigh  scattering  is  anisotropic. 
Once  we  get  deep  into  the  medium,  enough  scatterings  have  occurred  that 
photons  lose  memory  of  their  direction  of  entrance  and  thus  take  on  an 
isotropic  character. 

In  Table  2 we  considered  the  case  where  a dielectric  interface  of 
refractive  index  1.338  is  placed  at  tj  = 0.01  from  the  top  of  the  medium 
of  total  optical  thickness  t „ = 1000.  The  first  thing  to  note  is  that 
both  Hy  and  HQ  increase  by  roughly  a factor  of  two  across  the  interface 
with  the  exception  of  gQ  = 0.0199.  The  reason  for  this  anomaly  is  that 
the  direct  beam  is  the  major  contributer  to  the  downward  irradiance  just 
above  the  interface  and  at  this  low  sun  angle  most  of  the  radiation  is 
specularly  reflected  upwards  rather  than  being  transmitted  through  the 
interface.  For  the  higher  sun  angles  a substantial  part  of  the  direct 
radiation  is  transmitted  through  the  interface.  Another  factor  to  be 
considered  is  that  the  upwelling  radiation  beneath  the  interface,  which  lies 
outside  the  acceptance  cone,  gets  internally  reflected  and  therefore  the 
interface,  in  a sense,  acts  like  a trap  which  only  allows  a certain 
portion  of  the  upwelling  radiation  to  escape  while  the  remainder  gets 
continually  recycled.  Another  interesting  feature  of  this  case  is  that 
hy/Hy  at  t=0  and  just  above  the  interface  shows  significant  departure 
from  the  case  where  no  interface  is  present.  This  is  due  to  the  fact 
that  we  are  now  getting  a significant  specular  component  in  the  upwelling 
radiation  which  tends  to  drive  this  ratio  more  towards  1/SQ.  Also  by 
the  time  t=  10  is  reached  the  radiation  field  becomes  isotropic  once 
again;  however  the  irradiances  still  maintain  higher  values  at  equivalent 


optical  depths  compared  to  the  case  where  there  is  no  interface. 


In  Table  3 we  present  results  for  the  case  where  the  interface  is  moved 
to  tj  = 0.1.  For  this  case  Hjj  and  Hp  again  increase  across  the  interface 
for  all  values  of  z . The  anomaly  now  disappears  for  low  sun  angles  since 
the  major  contribution  to  the  downward  irradiance  just  above  the  interface 
is  due  to  the  diffuse  component  and  not  the  direct  beam,  resulting  in  more 
radiation  transmitted  through  the  interface.  It  should  also  be  noted  that 
ratios  of  quantities  are  basically  the  same  at  corresponding  optical  depths 
below  the  interface  regardless  of  their  position. 

In  Table  4 the  interface  has  been  moved  to  ij  = 10.0.  At  this  optical 
depth  the  effect  of  the  interface  can  not  be  seen  at  the  top  of  the  medium 
since  all  the  computed  quantities  are  essentially  the  same  compared  to  the 
case  where  no  interface  is  present.  Again  we  note  the  irradiance  increase 
across  the  interface  for  all  values  of  ZQ  and  then  the  monotonic  decrease 
as  we  move  to  larger  optical  depths. 

As  an  additional  aid  to  understanding  the  results  just  discussed,  we 
can  consider  Figs.  3-5.  Here  we  have  plotted  the  diffuse  component  of 
the  azimuthal ly  averaged  radiance  i.e.,  all  radiation  which  has  undergone 
at  least  one  scattering,  for  the  case  where  the  interface  has  been  placed 
at  ij  = 0.01.  In  Fig.  3a  we  show  the  upward  diffuse  component  of  radiance 
as  a function  of  m for  various  optical  depths.  The  results  have  been 
normalized  to  the  maximum  value  of  the  radiance  at  each  optical  depth.  We 
first  note  that  both  at  the  top  and  just  above  the  interface  the  upward 
diffuse  radiance  is  darkened  towards  the  horizon  but  as  we  cross  the 
interface  its  character  changes  markedly.  Now  as  we  move  deeper  into  the 
medium  where  single  scattering  effects  become  negligible  the  radiation  field 
becomes  isotropic.  In  Fig.  3b  we  show  the  diffuse  downward  radiance  for 
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the  same  £;  . Just  above  the  interface  the  radiance  is  highly  peaked  near 
the  horizon.  To  visualize  the  contribution  at  an  interior  point  one  must 
realize  that  for  the  downward  radiance  there  are  two  components  to  consider. 

First  there  is  transmission  of  the  external  source  through  the  layer  above 
the  detector  and  then  there  is  reflection  of  the  upwelling  radiance  from 
the  same  layer  above  the  detector.  Similar  reasoning  also  applies  to  the 
upward  radiance  at  an  interior  point.  Now  as  we  cross  the  interface  the 
largest  contribution  to  the  downward  radiance  is  coming  from  the  radiation 
traveling  upward  below  the  interface.  Therefore  for  0 < p < pc  the  upwelling 
radiance  is  totally  internally  reflected  in  the  downward  direction  whereas 
for  uc  < p < 1 the  upwelling  radiance  is  modified  by  the  reflection  coefficient  whicf 
is  close  to  unity  for  p - p but  falls  rapidly  for  p > p . Thus  the  rapid 
fall -off  for  p > pc  is  simply  due  to  the  reflectance  properties  of  the  inter- 
face. Again  we  notice  as  we  move  deeper  into  the  medium  this  strong  asymmetry 
gradually  disappears  and  when  we  reach  t = 10  the  radiation  field  is  isotropic. 

In  Figs.  4a  and  4b  we  give  similar  results  for  £ = 0.592.  The  results  are 
qualitatively  similar  to  the  previous  case.  Now  in  Figs.  5a  and  5b  we  have 
a very  low  sun  angle;  namely  = 0.0199  and  we  see  a radical  change  in  the 
distribution  of  the  radiation.  The  upward  diffuse  radiance  at  the  top  of 
the  atmosphere  is  now  a maximum  at  the  horizon.  This  phenomenon  is  due 
primarily  to  two  effects.  First  we  note  that  the  downward  radiance  just 
above  the  interface  (Fig.  5b)  is  peaked  near  the  horizon.  This  radiation  is 
specularly  reflected  into  the  upward  stream  since  the  reflection  coefficient 
is  high  for  grazing  incidence.  Secondly,  there  is  more  scattered  radiation 
reaching  the  detector  which  has  never  penetrated  the  interface  due  to  the 
larger  optical  path  lengths  along  the  directions  of  incidence  and  emergence. 

Also  for  the  detectors  just  beneath  the  interface  there  is  less  separation 
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in  the  downward  radiance  ratio  between  the  horizon  and  zenith  than  for  the 
cases  of  higher  sun  angles. 


CONCLUSION 

We  have  demonstrated  that  the  matrix  operator  technique  can  be  used  to 
give  highly  accurate  results  for  the  equation  of  transfer  when  coupling  an 
atmosphere  and  ocean  system.  The  advantages  of  the  method  are  as  follows: 

1.  Results  for  all  solar  incident  angles  are  obtained  in  a single 
computer  run. 

2.  Any  type  of  ocean  bottom  can  be  used  provided  one  knows  its 
reflection  properties.  If  one  assumes  a Lambert  type  reflecting 
surface  then  results  can  be  obtained  for  many  different  bottom 
albedos  in  a single  computer  run. 

3.  One  can  easily  incorporate  internal  detectors  where  all  relevant 
physically  measurable  quantities  can  be  calculated. 

4.  Large  numbers  of  inhomogenous  layers  can  be  coupled  to  give  a 
completely  realistic  representation  of  the  atmosphere-ocean  system. 

5.  One  can  also  easily  obtain  results,  from  the  same  computer  run,  for 
an  ocean  whose  depth  can  be  altered  by  successively  placing  the 
bottom  at  each  internal  detector  location,  therefore  one  can  explore 
bottom  effects  which  are  relevant  to  coastal  zones. 

6.  Computations  are  quite  fast.  For  each  case  presented  in  the  paper 
the  computation  time  was  approximately  30  sec  on  an  Amdahl  470V-6. 
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(Table  1.  Computed  optical  quantities  as  a function  of  optical  depth  x for  a homoge- 
nous system  of  optical  thickness  xmax  = 1000.  Rayleigh  scattering  is  assumed  every- 
where and  aiQ  = 1.0.  The  input  solar  flux  is  assumed  to  be  unity  in  a plane  perpen- 
dicular to  the  incoming  beam. 


T 

^0 

HU 

hd 

hd  ' HU 

hu/hd 

VHu 

VH0 

hU  + hD 

0 

0.980 

0.978 

0.980 

1.63-3 

0.998 

1.87 

1.02 

2.83 

0 

0.592 

0.591 

0.592 

7.42-4 

0.999 

2.00 

1.69 

2.18 

0 

0.0199 

1 . 98-2* 

1.99-2 

1.18-5 

0.999 

3.28 

50.4 

1.06 

.01 

0.980 

0.983 

0.984 

1.63-3 

0.998 

1.87 

1.07 

2.90 

.01 

0.592 

0.592 

0.592 

7.42-4 

0.999 

2.00 

1.75 

2.22 

.01 

0.0199 

1.61-2 

1.61-2 

1.18-5 

0.999 

3.00 

38.8 

0.675 

0.1 

0.980 

1.01 

1.01 

1.63-3 

0.998 

1.90 

1.30 

3.24 

0.1 

0.592 

0.595 

0.596 

7.42-4 

0.999 

2.01 

1.96 

2.36 

0.1 

0.0199 

9.96-3 

9.97-3 

1.18-5 

0.999 

2.10 

3.59 

5.67-2 

1.0 

0.980 

1.15 

1.16 

1.63-3 

0.999 

1.97 

1.84 

4.40 

1.0 

0.592 

0.576 

0.576 

7.42-4 

0.999 

2.02 

2.09 

2.37 

1.0 

0.0199 

8.97-3 

8.99-3 

1 .18-5 

0.999 

2.01 

2.06 

3.65-2 

10 

0.980 

1.21 

1.21 

1 .63-3 

0.999 

2.00 

2.00 

4.84 

10 

0.592 

0.551 

0.552 

7.42-4 

0.999 

2.00 

2.00 

2.21 

10 

0.0199 

8.76-3 

8.78-3 

1.18-5 

0.999 

2.00 

2.00 

3.51-2 

20 

0.980 

1.20 

1.20 

1.63-3 

0.999 

2.00 

2.00 

4.80 

20 

0.592 

0.546 

0.546 

7.42-4 

0.999 

2.00 

2.00 

2.18 

20 

0.0199 

8.68-3 

8.69-3 

1.18-5 

0.999 

2.00 

2.00 

3.47-2 

500 

0.980 

0.611 

0.613 

1.63-3 

0.997 

2.00 

2.00 

2.45 

500 

0.592 

0.278 

0.279 

7.42-4 

0.997 

2.00 

2.00 

1.11 

500 

0.0199 

4.43-3 

4.44-3 

1.18-5 

0.997 

2.00 

2.00 

1.77-2 

* * 

The  notation  -2  indicates  the  power  of  ten  to  be  multiplied  by  the  preceeding  number. 
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Table  2.  Computed  optical  quantities  as  a function  of  t for  a medium  of  total  optical 
thickness  xmax  = 1000  with  a dielectric  interface  at  ij  = 0.01.  Rayleigh  scattering 
is  assumed  both  above  and  below  the  interface  and  uq  = 1.0.  The  input  solar  flux  is 
assumed  to  be  unity  in  a plane  perpendicular  to  the  incoming  beam. 


T 

^0 

Hu 

hd 

VHU 

VHo 

hu/Hu 

hD/HD 

hU  + hD 

0 

0.980 

0.978 

0.980 

2.57-3 

0.997 

1.79 

1.02 

2.75 

0 

0.592 

0.590 

0.592 

1.42-3 

0.998 

1.83 

1 .69 

2.08 

0 

0.0199 

1.98-2 

1.99-2 

1.60-5 

0.999 

19.1 

50.4 

1.38 

0.01A* 

0.980 

0.981 

0.984 

2.57-3 

0.997 

1.78 

1 .07 

2.80 

0.01A 

0.592 

0.591 

0.592 

1.42-3 

0.998 

1.81 

1.74 

2.10 

0.01A 

0.0199 

1.81-2 

1.81-2 

1.60-5 

0.999 

31.0 

35.1 

1.20 

0.01B 

0.980 

1 .82 

1.82 

2.57-3 

0.999 

1 .98 

1.88 

7.03 

0.01B 

0.592 

1.09 

1.09 

1.42-3 

0.999 

2.01 

2.04 

4.41 

0.01B 

0.0199 

1.24-2 

1.24-2 

1.60-5 

0.999 

2.01 

2.08 

5.10-2 

0.1 

0.980 

1.83 

1.83 

2.57-3 

0.999 

1.98 

1.89 

7.07 

0.1 

0.592 

1.09 

1.09 

1.42-3 

0.999 

2.01 

2.04 

4.40 

0.1 

0.0199 

1.24-2 

1.24-2 

1.60-5 

0.999 

2.01 

2.08 

5.08-2 

1.0 

0.980 

1 .88 

1.88 

2.57-3 

0.999 

1.99 

1.94 

7.41 

1 .0 

0.592 

1.07 

1.08 

1 .42-3 

0.999 

2.00 

2.02 

4.33 

1.0 

0.0199 

1.21-2 

1 .21-2 

1.60-5 

0.999 

2.01 

2.03 

4.90-2 

10 

0.980 

1.91 

1.91 

2.57-3 

0.999 

2.00 

2.00 

7.64 

10 

0.592 

1.05 

1.06 

1.42-3 

0.999 

2.00 

2.00 

4.22 

10 

0.0199 

1.19-2 

1.19-2 

1.60-5 

0.999 

2.00 

2.00 

4.75-2 

20 

0.980 

1.89 

1.89 

2.57-3 

0.999 

2.00 

2.00 

7.56 

20 

0.592 

1.04 

1.04 

1.42-3 

0.999 

2.00 

2.00 

4.18 

20 

0.0199 

1.17-2 

1.18-2 

1.60-5 

0.999 

2.00 

2.00 

4.70-2 

500 

0.980 

0.964 

0.966 

2.57-3 

0.997 

2.00 

2.00 

3.86 

500 

0.592 

0.533 

0.534 

1.42-3 

0.997 

2.00 

2.00 

2.13 

500 

0.0199 

5.99-3 

6.01-3 

1.60-5 

0.997 

2.00 

2.00 

2.40-2 

^Superscripts  A and  11  refer  to  detectors  placed  just  above  and  just  below  the  interface 
respectively. 


Table 

3 . Same 

as  Table  2 except  Tj 

= 0.1. 

T 

Co 

Hu 

hd 

hd'  hu 

Vhd 

VHu 

*Vhd 

VhD 

0 

0.980 

0.977 

0.980 

2.60-3 

0.997 

1.85 

1.02 

2.81 

0 

0.592 

0.590 

0.592 

1.41-3 

0.997 

1.95 

1.69 

2.15 

0 

0.0199 

1.98-2 

1.99-2 

2.21-5 

0.999 

3.35 

50.4 

1.07 

0.01 

0.980 

0.981 

0.984 

2.60-3 

0.997 

1.86 

1.07 

2.88 

0.01 

0.592 

0.591 

0.592 

1.41-3 

0.998 

1.95 

1.75 

2.19 

0.01 

0.0199 

1.61-2 

1.61-2 

2.21-5 

0.999 

3.08 

3.88 

0.676 

0.1A 

0.980 

1.01 

1.01 

2.60-3 

0.997 

1.87 

1.29 

3.20 

0.1A 

0.592 

0.593 

0.594 

1.41-3 

0.998 

1 .93 

1.95 

2.30 

0.1A 

0.0199 

1.01-2 

1.01-2 

2.21-5 

0.998 

2.98 

3.60 

6.65-2 

0.1B 

0.980 

1.86 

1.86 

2.60-3 

0.999 

1 .98 

1 .90 

7.21 

0.1B 

0.592 

1.08 

1.08 

1.41-3 

0.999 

2.01 

2.04 

4.37 

0.1B 

0.0199 

1.69-2 

1.69-2 

2.21-5 

0.999 

2.01 

2.04 

6.84-2 

1.0 

0.980 

1.91 

1.91 

2.60-3 

0.999 

1.99 

1.95 

7.52 

1.0 

0.592 

1 .07 

1.07 

1.41-3 

0.999 

2.00 

2.02 

4.30 

1.0 

0.0199 

1.67-2 

1.67-2 

2.21-5 

0.999 

2.00 

2.02 

6.71-2 

10 

0.980 

1.93 

1.94 

2.60-3 

0.999 

2.00 

2.00 

7.74 

10 

0.592 

1.05 

1.05 

1.41-3 

0.999 

2.00 

2.00 

4.19 

10 

0.0199 

1.64-2 

1.64-2 

2.21-5 

0.999 

2.00 

2.00 

6.56-2 

20 

0.980 

1.91 

1.92 

2.60-3 

0.999 

2.00 

2.00 

7.66 

20 

0.592 

1.03 

1.04 

1.41-3 

0.999 

2.00 

2.00 

4.14 

20 

0.0199 

1.62-2 

1.62-2 

2.21-5 

0.999 

2.00 

2.00 

6.49-2 

500 

0.980 

0.976 

0.979 

2.60-3 

0.997 

2.00 

2.00 

3.91 

500 

0.592 

0.528 

0.530 

1.41-3 

0.997 

2.00 

2.00 

2.12 

500 

0.0199 

8.27-3 

8.30-3 

2.21-5 

0.997 

2.00 

2.00 

3.31-5 
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Table 

4 . Same 

as  Table  2 except  ij 

= 10.0. 

T 

^0 

Hu 

H0 

VHU 

hij/hd 

VHu 

hD/HD 

hU  + hD 

0 

0.980 

0.977 

0.980 

2.89-3 

0.997 

1.87 

1.02 

2.83 

0 

0.592 

0.590 

0.592 

1.32-3 

0.998 

2.00 

1.69 

2.18 

0 

0.0199 

1.98-2 

1.99-2 

2.09-5 

0.999 

3.28 

50.4 

1.06 

0.01 

0.980 

0.981 

0.984 

2.89-3 

0.997 

1.87 

1.07 

2.89 

0.01 

0.592 

0.591 

0.592 

1.32-3 

0.998 

2.00 

1.75 

2.22 

0.01 

0.0199 

1.61-2 

1.61-2 

2.09-5 

0.999 

3.00 

38.8 

0.675 

0.1 

0.980 

1.01 

1.01 

2.89-3 

0.997 

1.90 

1.30 

3.24 

0.1 

0.592 

0.595 

0.596 

1.32-3 

0.998 

2.01 

1.96 

2.36 

0.1 

0.0199 

9.45-3 

9.97-3 

2.09-5 

0.998 

2.10 

3.59 

5.67-2 

1.0 

0.980 

1.15 

1.15 

2.89-3 

0.997 

1.97 

1.84 

4.39 

1.0 

0.592 

0.575 

0.576 

1.32-3 

0.998 

2.02 

2.09 

2.36 

1.0 

0.0199 

8.96-3 

8.98-3 

2.09-5 

0.998 

2.01 

2.06 

3.65-2 

10A 

0.980 

1.20 

1.20 

2.89-3 

0.998 

2.00 

2.00 

4.80 

10A 

0.592 

0.546 

0.548 

1.32-3 

0.998 

2.00 

2.00 

2.19 

10A 

0.0199 

8.69-3 

8.71-3 

2.09-5 

0.998 

2.00 

2.00 

3.48-2 

10B 

0.980 

2.15 

2.15 

2.89-3 

0.999 

2.00 

2.00 

8.59 

10B 

0.592 

0.978 

0.980 

1.32-3 

0.999 

2.00 

2.00 

3.91 

10B 

0.0199 

1.55-2 

1.56-2 

2.09-5 

0.999 

2.00 

2.00 

6.22-2 

20 

0.980 

2.13 

2.13 

2.89-3 

0.999 

2.00 

2.00 

8.51 

20 

0.592 

0.968 

0.970 

1.32-3 

0.999 

2.00 

2.00 

3.88 

20 

0.0199 

1.54-2 

1.54-2 
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FIGURE  CAPTIONS 

Fig.  1.  Geometric  representation  of  angular  variables  used  in  the  atmosphere- 
ocean  system. 

Fig.  2.  a)  Function  transformation  used  for  mapping  ocean  angular  variables 
into  atmospheric  angular  variables,  b)  Function  transformation  used 
to  map  atmospheric  angular  variables  into  ocean  angular  variables. 

Fig.  3.  a)  Upward  diffuse  azimuthal ly  averaged  radiance  1°  versus  y for  a 
conservative  Rayleigh  scattering  medium  of  total  optical  thickness 

= 1000  with  a dielectric  interface  at  t,  = 0.01  for  optical 
max  I 

depth  values  of  0,  0.01A,  0.01B,  1.0,  and  10.0.  Superscripts  A and 
B refer  to  detectors  just  above  and  just  below  the  interface  res- 
pectively. The  solar  incident  angle  is  £ = 0.980  (a  = 1.5°)  and 
the  results  are  normalized  to  their  maximum  value  at  each  optical 
depth,  b)  Same  as  a)  except  results  are  for  the  downward  diffuse  1°. 

Fig.  4.  a)  Same  as  3a  except  £ = 0.592  (<»  = 53.7°).  b)  Same  as  3b  except 
f = 0.592  (u  = 53.7°). 

a)  Same  as  3a  except  = 0.0199  (a  = 88.9°).  b)  Same  as  3b  except 
£0  = 0.0199  (a  = 88.9°). 


Fig.  5. 
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